Liquid argon benchmarks

Sebastian Micluța-Câmpeanu, Mikhail Vaganov

The purpose of these benchmarks is to compare several integrators for use in molecular dynamics simulation. We will use a simulation of liquid argon form the examples of NBodySimulator as test case.

using ProgressLogging
using NBodySimulator, OrdinaryDiffEq, StaticArrays
using Plots, DataFrames, StatsPlots

function setup(t)
    T = 120.0 # K
    kb = 1.38e-23 # J/K
    ϵ = T * kb # J
    σ = 3.4e-10 # m
    ρ = 1374 # kg/m^3
    m = 39.95 * 1.6747 * 1e-27 # kg
    N = 216
    L = (m*N/ρ)^(1/3)
    R = 2.25σ
    v_dev = sqrt(kb * T / m) # m/s

    _L = L / σ
     = 1.0
     = 1.0
    _m = 1.0
    _v = v_dev / sqrt(ϵ / m)
    _R = R / σ

    bodies = generate_bodies_in_cell_nodes(N, _m, _v, _L)
    lj_parameters = LennardJonesParameters(, , _R)
    pbc = CubicPeriodicBoundaryConditions(_L)
    lj_system = PotentialNBodySystem(bodies, Dict(:lennard_jones => lj_parameters));
    simulation = NBodySimulation(lj_system, (0.0, t), pbc, /T)

    return simulation
end
setup (generic function with 1 method)

In order to compare different integrating methods we will consider a fixed simulation time and change the timestep (or tolerances in the case of adaptive methods).

function benchmark(energyerr, rts, bytes, allocs, nt, nf, t, configs)
    simulation = setup(t)
    prob = SecondOrderODEProblem(simulation)
    for config in configs
        alg = config.alg
        sol, rt, b, gc, memalloc = @timed solve(prob, alg();
            save_everystep=false, progress=true, progress_name="$alg", config...)
        result = NBodySimulator.SimulationResult(sol, simulation)
        ΔE = total_energy(result, t) - total_energy(result, 0)
        energyerr[alg] = ΔE
        rts[alg] = rt
        bytes[alg] = b
        allocs[alg] = memalloc
        nt[alg] = sol.destats.naccept
        nf[alg] = sol.destats.nf + sol.destats.nf2
    end
end

function run_benchmark!(results, t, integrators, tol...; c=ones(length(integrators)))
    @progress "Benchmark at t=$t" for τ in zip(tol...)
        runtime = Dict()
        ΔE = Dict()
        nt = Dict()
        nf = Dict()
        b = Dict()
        allocs = Dict()
        cfg = config(integrators, c, τ...)

        GC.gc()
        benchmark(ΔE, runtime, b, allocs, nt, nf, t, cfg)
        get_tol(idx) = haskey(cfg[idx], :dt) ? cfg[idx].dt : (cfg[idx].abstol, cfg[idx].rtol)

        for (idx,i) in enumerate(integrators)
            push!(results, [string(i), runtime[i], get_tol(idx)..., abs(ΔE[i]), nt[i], nf[i], c[idx]])
        end
    end
    return results
end
run_benchmark! (generic function with 1 method)

We will consider symplectic integrators first

symplectic_integrators = [
    VelocityVerlet,
    #VerletLeapfrog,
    PseudoVerletLeapfrog,
    McAte2,
    #CalvoSanz4,
    #McAte5,
    Yoshida6,
    #KahanLi8,
    #SofSpa10
]
4-element Array{DataType,1}:
 VelocityVerlet
 PseudoVerletLeapfrog
 McAte2
 Yoshida6

Let us run a short simulation to see the cost per timestep for each method

config(integrators, c, τ) = [ (alg=a, dt=τ*cₐ) for (a,cₐ) in zip(integrators, c)]

t = 35.0
τs = 1e-3

# warmup
c_symplectic = ones(length(symplectic_integrators))
benchmark(Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), 10.,
    config(symplectic_integrators, c_symplectic, τs))

results = DataFrame(:integrator=>String[], :runtime=>Float64[], =>Float64[],
    :EnergyError=>Float64[], :timesteps=>Int[], :f_evals=>Int[], :cost=>Float64[]);
run_benchmark!(results, t, symplectic_integrators, τs)

4 rows × 7 columns

integratorruntimeτEnergyErrortimestepsf_evalscost
StringFloat64Float64Float64Int64Int64Float64
1VelocityVerlet58.16890.0010.0036656235000700021.0
2PseudoVerletLeapfrog58.75720.0010.0208884350001050021.0
3McAte257.13420.0010.00980374350001050021.0
4Yoshida6227.8940.0010.0219439350005250021.0

The cost of a timestep can be computed as follows

c_symplectic .= results[!, :runtime] ./ results[!, :timesteps]
c_Verlet = c_symplectic[1]
c_symplectic /= c_Verlet
4-element Array{Float64,1}:
 1.0
 1.010114030346147
 0.982211928654614
 3.9178041722138084

were we have normalized the cost to the cost of a VelocityVerlet step.

Let us now benchmark the solvers for a fixed simulation time and variable timestep

t = 45.0
τs = 10 .^range(-4, -3, length=10)

results = DataFrame(:integrator=>String[], :runtime=>Float64[], =>Float64[],
    :EnergyError=>Float64[], :timesteps=>Int[], :f_evals=>Int[], :cost=>Float64[]);
run_benchmark!(results, t, symplectic_integrators, τs, c=c_symplectic)

40 rows × 7 columns

integratorruntimeτEnergyErrortimestepsf_evalscost
StringFloat64Float64Float64Int64Int64Float64
1VelocityVerlet771.3010.00010.002053744500009000021.0
2PseudoVerletLeapfrog758.320.0001010110.0021517144549513364871.01011
3McAte2773.2769.82212e-50.00027167945815013744520.982212
4Yoshida6786.3250.000391780.0034387211486117229173.9178
5VelocityVerlet603.3110.0001291550.003468873484196968401.0
6PseudoVerletLeapfrog591.6410.0001304610.0034931334493110347951.01011
7McAte2610.0590.0001268580.0033625335472910641890.982212
8Yoshida6611.4170.0005060040.001760998893313339973.9178
9VelocityVerlet468.3430.000166810.002438572697685395381.0
10PseudoVerletLeapfrog452.9270.0001684970.007898532670678012031.01011
11McAte2463.0710.0001638430.00153742746548239640.982212
12Yoshida6470.5270.0006535290.002974776885710328573.9178
13VelocityVerlet361.3190.0002154430.01085092088724177461.0
14PseudoVerletLeapfrog356.2870.0002176220.004854932067816203451.01011
15McAte2361.2060.0002116110.006134042126556379670.982212
16Yoshida6359.0160.0008440650.0677113533147997123.9178
17VelocityVerlet276.6630.0002782560.004961421617223234461.0
18PseudoVerletLeapfrog276.2590.000281070.002077051601034803111.01011
19McAte2285.0250.0002733060.0006883711646514939550.982212
20Yoshida6284.420.001090150.0185529412796191873.9178
21VelocityVerlet213.5360.0003593810.0003693991252162504341.0
22PseudoVerletLeapfrog211.6680.0003630160.005589551239623718881.01011
23McAte2218.8720.0003529890.005232211274833824510.982212
24Yoshida6218.5180.001407990.00593719319614794173.9178
25VelocityVerlet164.2440.0004641590.00704297969501939021.0
26PseudoVerletLeapfrog163.670.0004688530.00314934959792879391.01011
27McAte2171.9250.0004559020.00749259987062961200.982212
28Yoshida6169.2970.001818480.000755624247463711923.9178
29VelocityVerlet129.6790.0005994840.0267297750651501321.0
30PseudoVerletLeapfrog127.1690.0006055470.00661891743132229411.01011
31McAte2130.7230.0005888210.0105314764242292740.982212
32Yoshida6130.9370.002348660.0323632191602874023.9178
33VelocityVerlet100.9030.0007742640.0220466581201162421.0
34PseudoVerletLeapfrog98.39570.0007820950.0263143575381726161.01011
35McAte2100.6970.0007604910.00945157591731775210.982212
36Yoshida6103.340.003033410.033734148352225273.9178
37VelocityVerlet77.30160.0010.012826845001900041.0
38PseudoVerletLeapfrog75.67410.001010110.00695002445501336521.01011
39McAte278.54040.0009822120.00132507458151374470.982212
40Yoshida680.92940.00391780.0533705114871723073.9178

The energy error as a function of runtime is given by

@df results plot(:EnergyError, :runtime, group=:integrator,
    xscale=:log10, yscale=:log10, xlabel="Energy error", ylabel="Runtime (s)")

Looking at the runtime as a function of timesteps, we can observe that we have a linear dependency for each method, and the slope is the previously computed cost per step.

@df results plot(:timesteps, :runtime, group=:integrator,
    xscale=:log10, yscale=:log10, xlabel="Number of timesteps", ylabel="Runtime (s)")

We can also consider a longer simulation time

t = 100.0

τs = 10 .^range(-4, -3, length=5)

results = DataFrame(:integrator=>String[], :runtime=>Float64[], =>Float64[],
    :EnergyError=>Float64[], :timesteps=>Int[], :f_evals=>Int[], :cost=>Float64[]);
#run_benchmark!(results, t, symplectic_integrators, τs, c=c_symplectic)

The energy error as a function of runtime is given by

#@df results plot(:EnergyError, :runtime, group=:integrator,
#    xscale=:log10, yscale=:log10, xlabel="Energy error", ylabel="Runtime (s)")

We can also look at the energy error history

function benchmark(energyerr, rts, ts, t, configs)
    simulation = setup(t)
    prob = SecondOrderODEProblem(simulation)
    for config in configs
        alg = config.alg
        sol, rt = @timed solve(prob, alg(); progress=true, progress_name="$alg", config...)
        result = NBodySimulator.SimulationResult(sol, simulation)
        ΔE(t) = total_energy(result, t) - total_energy(result, 0)
        energyerr[alg] = [ΔE(t) for t in sol.t[2:end]]
        rts[alg] = rt
        ts[alg] = sol.t[2:end]
    end
end

ΔE = Dict()
rt = Dict()
ts = Dict()
configs = config(symplectic_integrators, c_symplectic, 2.3e-4)
benchmark(ΔE, rt, ts, 45., configs)

plt = plot(xlabel="Rescaled Time", ylabel="Energy error", legend=:bottomleft);
for c in configs
    plot!(plt, ts[c.alg], abs.(ΔE[c.alg]), label="$(c.alg), $(rt[c.alg])s", yscale=:log10)
end
plt

Now, let us compare some adaptive methods

adaptive_integrators=[
    # DPRKN
    DPRKN6,
    DPRKN8,
    DPRKN12,
    # others
    Tsit5,
    Vern7,
    Vern9
]

config(integrators, c, at, rt) = [ (alg=a, abstol=at, rtol=rt) for a in integrators]

t = 45.0
ats = 10 .^range(-20, -14, length=10)
rts = 10 .^range(-20, -14, length=10)

# warmup
benchmark(Dict(), Dict(), Dict(), Dict(), Dict(), Dict(), 10.,
    config(adaptive_integrators, 1, ats[1], rts[1]))

results = DataFrame(:integrator=>String[], :runtime=>Float64[], :abstol=>Float64[],
    :reltol=>Float64[], :EnergyError=>Float64[], :timesteps=>Int[], :f_evals=>Int[], :cost=>Float64[]);
run_benchmark!(results, t, adaptive_integrators, ats, rts)

60 rows × 8 columns

integratorruntimeabstolreltolEnergyErrortimestepsf_evalscost
StringFloat64Float64Float64Float64Int64Int64Float64
1DPRKN642.5541.0e-201.0e-200.04421657661571871.0
2DPRKN895.59011.0e-201.0e-200.0566201118371248941.0
3DPRKN12187.8471.0e-201.0e-200.00742181125622295761.0
4Tsit524.24491.0e-201.0e-200.03337574234281131.0
5Vern746.11591.0e-201.0e-200.08628395123532621.0
6Vern9212.8231.0e-201.0e-200.107903143772460821.0
7DPRKN643.45524.64159e-204.64159e-200.01116997778585661.0
8DPRKN895.94764.64159e-204.64159e-200.00555696117381240341.0
9DPRKN12185.84.64159e-204.64159e-200.0348998125812299001.0
10Tsit524.51814.64159e-204.64159e-200.2491014206280111.0
11Vern746.69124.64159e-204.64159e-200.01922285147536421.0
12Vern9214.0794.64159e-204.64159e-200.151188142612442901.0
13DPRKN641.882.15443e-192.15443e-190.05611547604567531.0
14DPRKN894.65992.15443e-192.15443e-190.00726108117181232241.0
15DPRKN12185.6952.15443e-192.15443e-190.00245804125652294681.0
16Tsit524.90252.15443e-192.15443e-190.8028254256286411.0
17Vern746.97522.15443e-192.15443e-190.02008285177540321.0
18Vern9212.7112.15443e-192.15443e-190.132925142502443221.0
19DPRKN642.74881.0e-181.0e-180.03194397676576141.0
20DPRKN895.41991.0e-181.0e-180.0344401117551239441.0
21DPRKN12186.9581.0e-181.0e-180.00742968126122302061.0
22Tsit524.09271.0e-181.0e-181.322274190279091.0
23Vern746.58461.0e-181.0e-180.02663895121534221.0
24Vern9219.3111.0e-181.0e-180.138736145832515221.0
25DPRKN642.35174.64159e-184.64159e-180.03206987711578241.0
26DPRKN893.40784.64159e-184.64159e-180.0369092116301216641.0
27DPRKN12188.3924.64159e-184.64159e-180.0420185126202307101.0
28Tsit524.19524.64159e-184.64159e-181.25364176278851.0
29Vern745.59414.64159e-184.64159e-180.1148675083527321.0
30Vern9212.9724.64159e-184.64159e-180.149868142792451061.0
31DPRKN641.83942.15443e-172.15443e-170.06469667600567741.0
32DPRKN897.22442.15443e-172.15443e-170.00196891119891275541.0
33DPRKN12186.5782.15443e-172.15443e-170.049486126072303321.0
34Tsit524.3562.15443e-172.15443e-170.8092224205280771.0
35Vern744.86332.15443e-172.15443e-170.004877584990517621.0
36Vern9214.4152.15443e-172.15443e-170.156428142552448181.0
37DPRKN640.61891.0e-161.0e-160.009810897485554791.0
38DPRKN896.47971.0e-161.0e-160.0218561118361250541.0
39DPRKN12186.4441.0e-161.0e-160.0159265125302286221.0
40Tsit523.86951.0e-161.0e-160.544874222280831.0
41Vern745.1111.0e-161.0e-160.01574255075527621.0
42Vern9212.741.0e-161.0e-160.126142672445461.0
43DPRKN643.02534.64159e-164.64159e-160.06572887604569631.0
44DPRKN894.38794.64159e-164.64159e-160.0256069116051217541.0
45DPRKN12191.3064.64159e-164.64159e-160.0330368128052347781.0
46Tsit524.17734.64159e-164.64159e-160.0416914192277951.0
47Vern746.28814.64159e-164.64159e-160.03798665071529721.0
48Vern9221.1034.64159e-164.64159e-160.164754144002482261.0
49DPRKN642.46922.15443e-152.15443e-150.005848867646570331.0
50DPRKN895.41652.15443e-152.15443e-150.0241828117611237941.0
51DPRKN12185.4842.15443e-152.15443e-150.00616829125852297021.0
52Tsit524.29342.15443e-152.15443e-150.3688564179277171.0
53Vern746.30082.15443e-152.15443e-150.09363715147536621.0
54Vern9215.4232.15443e-152.15443e-150.169246143662469621.0
55DPRKN641.61951.0e-141.0e-140.09260447515559271.0
56DPRKN894.23631.0e-141.0e-140.029884116511223941.0
57DPRKN12187.8251.0e-141.0e-140.0207755126572311961.0
58Tsit524.39421.0e-141.0e-140.3569084192280111.0
59Vern745.66691.0e-141.0e-140.06985184993517021.0
60Vern9216.3311.0e-141.0e-140.153606142522445461.0

The energy error as a function of runtime is given by

@df results plot(:EnergyError, :runtime, group=:integrator,
    xscale=:log10, yscale=:log10, xlabel="Energy error", ylabel="Runtime (s)")

If we consider the number of function evaluations instead, we obtain

@df results plot(:EnergyError, :f_evals, group=:integrator,
    xscale=:log10, yscale=:log10, xlabel="Energy error", ylabel="Number of f evals")

We can also consider a longer simulation time

t = 100.0

ats = 10 .^range(-20, -14, length=10)
rts = 10 .^range(-20, -14, length=10)

results = DataFrame(:integrator=>String[], :runtime=>Float64[], :abstol=>Float64[],
    :reltol=>Float64[], :EnergyError=>Float64[], :timesteps=>Int[], :f_evals=>Int[], :cost=>Float64[]);
#run_benchmark!(results, t, integrators, ats, rts)

The energy error as a function of runtime is given by

#@df results plot(:EnergyError, :runtime, group=:integrator,
#    xscale=:log10, yscale=:log10, xlabel="Energy error", ylabel="Runtime (s)")